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I. INTRODUCTION 



During the last decade, much effort has been devoted to experimentally detecting the quark-gluon plasma state 
in high energy heavy-ion collisions. In order to extract an unambiguous signal of quark-gluon plasma production 
from heavy-ion collision experiments, theoretical understanding on the nature of the finite-temperature chiral phase 
transition and the thermodynamic properties of quark-gluon plasma is indispensable. In particular, the equation of 
state (EOS) belongs to the most basic category of information needed in phenomenological investigations of heavy-ion 
collisions. 

Extensive numerical studies have been pursued in lattice QCD to derive the equation of state from first principles [0 . 
Within the quenched approximation in which effects of dynamical quark-pair creation and annihilation are neglected, 
precise results have been established. Continuum extrapolations of the lattice results have been made with various 
lattice actions, finding a good agreement within errors of a few percent [||-||. For the pressure, a detailed comparison 
of the results from the integral method and the derivative method have also been made . The problem of non-zero 
pressure gap at the transition point with the derivative method has been solved by a non-perturbative calculation of 
anisotropy coefficients Q . 

The next step towards a realistic quark-gluon plasma simulation is to include dynamical quarks, which clearly play 
a significant role in the real world through chiral symmetry. Until recently, EOS with dynamical quarks has been 
computed only with the Kogut-Susskind (staggered) quark action or its improved form |p|-p^ . Strictly speaking, the 
staggered quark action only allows the number of fiavors to be a multiple of four, and providing a mass difference 
within the four-fold multiplet is not straightforward. It has also been found that the critical scaling for two-flavor QCD 
extracted with this formalism ||Tl|-[l^ does not agree with the theoretically expected 0(4) values. These features of the 
Kogut-Susskind quark action make it imperative that the quark-gluon plasma properties be explored with alternative 
quark actions. In this article we present results on the equation of state obtained with the Wilson quark action in an 
improved form. 

Study of finite temperature QCD with Wilson-type quark action has been difficult for two reasons. First, explicit 
chiral symmetry breaking complicates the phase diagram analysis []T5|-pT[, which is basic for obtaining the equation 
of state. In this connection, an important role played by the parity-flavor broken phase has been realized, and 
the phase structure for flnite temporal lattice sizes has been understood ||l9|-plt. 

Another difficulty has been that, when the standard plaquette gauge action and the standard Wilson quark action 
are adopted, the system exhibits severe lattice artifacts on coarse lattices with the temporal lattice size iV* = 4 and 6. 
For example, the finite-temperature transition strengthens at intermediate quark masses to a first-order transition for 
Nt — 6 [EgUT^, while it should weaken as quarks become heavier. In this regard, it has been shown that improvement 
of the gauge action is effective in reducing the lattice artifacts in finite temperature QCD. Furthermore, the critical 
scaling around the two-flavor chiral transition obtained for a renormalization-group (RG) improved gauge action is 
consistent with the expected 0(4) universality class at Nt = A [p2| . 

These advances indicate that thermodynamic studies with Wilson-type quark actions are feasible if improved actions 
are employed. We have thus attempted a first calculation of EOS in QCD with two fiavors of dynamical quarks, using 
an RG improved gauge action coupled with a clover-improved Wilson quark action [Q. This combination of 
actions is motivated from our previous comparative study psf , where we found lattice discretization errors to be small 
with this action combination in both gluonic and hadronic observables at zero temperature. 

The phase structure and the critical temperature for this action combination have been studied in Ref. [p6t employing 
an iV( = 4 lattice. In this paper, we extend the study to an A't = 6 lattice. We then calculate the pressure and energy 
density for Nt = 4, and 6 employing the integral method We obtain EOS as a function of temperature for each 
fixed value of the renormalized quark mass, i.e., on each line of constant physics. We identify these lines by the ratio 
of pseudo-scalar and vector meson masses, mps/my, at zero temperature. Our results covers EOS over the range 
mps/wv — 0.65-0.95. 

The organiz atio n of the paper is as follows. Our lattice action and the simulation parameters are summarized in 



Sec. UJL In Sec. ID, we discuss the phase structure of QCD for our improved Wilson quark action at A''^ = 4 and 6. In 
Sec. |[V| the temperature scale and the lines of constant physics are studied. The RG beta functions, required in our 
calculation of EOS, are determined in Sec. 0. Results for EOS at A^t = 4 and 6 are presented in Sec. VI. Conclusions 



and discussions are given in Sec. VII 



II. ACTION AND SIMULATION PARAMETERS 

The gluon |23] and quark [e3 action we employ is defined by 
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S = Sg+ Sg, (1) 

= -Pi E col^^l^^(x) + E ^iK^^'i^)}^ (2) 



X. fi>iy 



f=l,2 x,y 

Here (3 = 6/5^, ci = —0.331, cq = 1 — 8ci and 

^'ic.y = Sxy - K^^Kl - '^f^)Ux,fiSa:+fi,y + (1 + 7t2)U^^^Sx,y+il} - SxyCswK^^ Cr^i/^)tiy (4) 



with the lattice field strength, 



where f^^ is the standard clover-shaped combination of gauge links. For the clover coefficient csw^ we adopt a 
mean-field value csw — P~'^^^ substituting the one-loop result for the plaquette P = 1 — 0.8412/3^^ |^^, which agrees 
within 8% with the values measured in our runs We also note that the one-loop result csw = 1 + 0.678//3 
p9[ is close to our choice csw = 1 + 0.631//3 -I- • • •. 

Our two-flavor simulation employs the standard HMC algorithm. Details of the algorithm implementation are the 
same as in Refs. [ p5|p^j2^ . The length of one trajectory is unity, and the molecular dynamics time step is chosen to 
yield an acceptance rate greater than about 80%. The inversion of quark matrix is made with the BiCGStab method. 
We measure the Polyakov loop and its susceptibility at every trajectory. Jack-knife errors of these expectation values 
are estimated with a bin size of 20-50 trajectories. Hadron propagators are measured at every 5 trajectories using 
point and exponentially smeared quark sources. 

In Ref. 12^], we studied the phase structure for our action combination on 16'^ x 4 lattices with a temporal lattice 
size Nt = 4. The simulation parameters are summarized in Table |[ The values of the hopping parameter K cover the 
range mps/wv ~ 0.60-0.98. In these simulations, we have also measured the observables required for EOS. We have 
since extended the simulation to an iVt = 6 lattice. Simulation parameters for Nt = 6 are compiled in Table H. For 
the spatial lattice size, we choose Ns ~ 16 both for Nt — 4: and 6. This enables us to commonly apply results obtained 
on a 16^ lattice to carry out the zero-temperature subtraction in the calculations of EOS, and to determine the lines 
of constant physics. For various tests, we also perform simulations on 8^ x 4 lattices as summarized in Appendix 



Parameters for the zero temperature simulations are compiled in Table III. On the zero temperature lattice 16 , we 
determine meson masses by a combined fit using both point and smeared sources assuming a double hyperbolic cosine 
form for the propagator. This procedure is necessitated by the fact that a plateau of effective mass is sometimes not 
quite clear due to a small temporal size of 16. Results for masses are summarized in Table |^ and plotted in Figs. |l| 
and |. 

In order to check the accuracy of mass results, we compare the values of mpg and my with those obtained in our 
previous high statistic simulations We find that, when Nta is smaller than about 6/mps {K — 0.1365 at /3 = 2.2 
and K — 0.1355 and 0.1360 at (3 = 2.3), several masses are inconsistent with those obtained on a 24^ x 48 lattice. In 
these cases, the effective mass on the 16'' lattice show no clear plateau up to the largest time separation 8, suggesting 
that the temporal lattice size of 16 is not sufficiently large to remove contamination from excited states. Therefore, 
we do not use data of mps and my on the 16** lattice when Nta < 6.0/mps. We instead include results from our 
previous study shown by open symbols in Figs. and|2[ in the analyses in the present study. 

Since the aspect ratio Ng/Nt = 16/6 = 2.666 • • • for iVt = 6 is smaller than Ng/Nt = 4 for Nt = 4, we also check 
the influence of the spatial volume on EOS. In the ideal gas limit of /? — > oo, analytic calculations as described in 
Appendix ^ show a 10% flnite size correction for Ns/Nt = 3 when Nt ~ 4-6 as compared to a 5% effect for Ng/Nt = 4. 
Perturbative estimates are not reliable close to the critical temperature, however. To study the spatial volume effect 
in this region, we make additional simulations at Ns/Nt = 2 on 8^ x 4 lattices and compare the results with those at 
Ns/Nt = 4 obtained on 16^ x 4 lattices. Details are presented in Appendix^. We flnd that the pressure at Ng/Nt = 2 
and 4 are consistent with each other within 1-2% except very near the critical temperature. We therefore conclude 
that finite volume corrections are reasonably controlled for our Nt = 6 lattices over the range of temperature we 
study. 
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III. PHASE STRUCTURE AND PSEUDO-CRITICAL TEMPERATURE 



Figure ^ summarize our results for the phase diagram. The solid line threading through filled circles denoted 
Kc{T = 0) is the location of the critical line where pion mass vanishes at zero temperature. It is the line of constant 
physics for massless quarks. Above the Kc{T = 0) line, parity-flavor symmetry of the Wilson-type quark action 
is broken spontaneously |l^j2^. Pion mass vanishes along the line since the pion becomes the massless mode of 
a second-order transition associated with this spontaneous breakdown. At zero temperature, the boundary of the 
parity- flavor broken phase (the Kc(T = 0) line) is expected to form a sharp cusp touching the free massless fermion 
point i^T = 1/8 at /3 = oo. 

For finite temporal sizes Nt corresponding to finite temperatures T = (aNt)^^, the parity-flavor broken phase 
retracts from the large P limit, forming a cusp ||l^,^. The boundary of the parity-flavor broken phase at Nt — 4, 
the Kc{Nt — 4) line, is shown by a thin line threading through open circles in Fig. ^. 

The dashed line Kt{Nt = 4) through open diamonds represents the finite-temperature pseudo-critical line for a 
temporal size Nt = 4, which is determined from the Polyakov loop and its susceptibility |^^. The region to the right 
of Kt (larger /3) is the high temperature quark-gluon plasma phase, and that to the left of Kt (smaller /3) is the 
low temperature hadron phase. The crossing point of the Kc{T = 0) and Kt{Nt = 4) lines is the finite-temperature 
chiral phase transition point ||l^ . As one observes in Fig. ^, the chiral transition point is located close to the cusp of 
the parity-flavor broken phase, with the difference expected to be 0{a). This is consistent with the picture that the 
massless pion, interpreted as the Nambu-Goldstone boson associated with spontaneous chiral symmetry breaking in 
the continuum limit, appears only in the cold phase. 

In Figs. ^ and |^ we present the expectation value of the Polyakov loop and its susceptibility obtained on a 16^ x 6 
lattice. We find a clear peak of the Polyakov loop susceptibility. Fitting the peak by a gaussian form using 3 or 4 
points around the peak, we determine Kt{Nt = 6) aX [3 — 2.0, 2.1, 2.2 and 2.3, as summarized in Table ^ The results 
are shown by filled diamonds denoted as Kt{Nt = 6) in Fig. |[ 

The pseudo-critical temperature Tpc in units of the zero-temperature vector meson mass was studied in Ref. p^ ] 
for Nt = 4. We repeat the study using the Nt = Q data. For this purpose, we interpolate the zero-temperature meson 
mass data to the Kt{Nt = 6) line by 
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(6) 

mya = By[-- — ]+Cy[-- — ] . (7) 



The values of mps/my and Tpc/my at Kt are summarized in Table | 

Results of Tpc/rriY as a function of (mps/my)^ are plotted in Fig. |6| for Nt — 4: and 6. We find that, in the 
range mp^/my — 0.65-0.95 we study, values of Tpc/my aX Nt — 4: and 6 agree within about 10%. 



IV. LINES OF CONSTANT PHYSICS 



In previous studies of EOS with staggered-type quark actions, the pressure and energy density are determined as 
functions of temperature for a fixed value of bare quark mass m'f^a. While rn'q^ a and Nt are practically easy to set 
in simulations, the system at different temperatures (or values of (3) will have different physical quark masses. This is 
not useful for phenomenological applications, and we need to evaluate the temperature dependence of thermodynamic 
observables for a fixed renormalized quark mass, i.e., on a line of constant physics. 

We identify the lines of constant physics in the parameter space (/3, K) by the values of mps/my at zero temperature. 
Deferring details of the interpolation procedure of hadron mass data to Sec. we show the lines of constant physics 
corresponding to the values mps/my = 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95 and 0.975 by solid lines in Fig. ^ In 
this figure, the bold line with open circles represents the critical line Kc{T = 0), corresponding to mpg/my = 0. The 
bold lines with open diamonds and triangles are the Kt lines for Nt = 4: and 6. 

We also attempt to determine the lines of constant temperature. Here, we adopt the pseudo-critical temperature Tpc 
on the same line of constant physics as the unit of temperature, where the temperature itself is determined through 
the zero-temperature vector meson mass mya as 

T f 

— iP, - X mva(/3, K) ' ^^"^ 
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The ratio Tpc/my for the pseudo-critical temperature Tpc is obtained by tuning 13 and K along the Kt line for each 
Nt- We then interpolate Tpc/my as a function of (mps/mv)^ by a Pade-type ansatz, 

7^ ^ ^ 1 + B (mps/TOy)^ 
my 1 + C (mpg/my)2 ■ 

We obtain A = 0.2253(71), B = -0.933(17), C = -0.820(39) with x^/Ndf = 1.61/5 for Nt = 4, and A = 0.261(19), 
B = -0.873(35), C = -0.624(108) with x^/Ndf = 0.74/1 for Nt = 6. Dashed and dot-dashed lines shown in Fig. | 
are the fit results for Nt = ^ and 6. 

The ansatz (||) does not incorporate the 0(4) scahng behavior rpc(wps) — Tpc(O) -I- cm1/g^ with 1/pS = 0.54 
expected close to the chiral limit or the constraint Tpc = in the heavy quark limit mps/my — 1. Fits satisfying 
these constraints may be attempted, e.g., by replacing (mps/my)^ with (m.pg/my)^/'''' and giving an additional 
factor (1 — (mps/my)^). They yield curves which are close to that of (^ but agree less well with data. Since we 
employ such fits for the purpose of interpolating the data for Tpc/my over the quark mass range mps/my ~ 0.65-0.95 
of our study, we choose to adopt (^) in the following analyses. 

For Nt = 6, since the data for Tpc/my covers only the range mps/iTiY = 0.725-0.972, we have to extrapolate the 
fit result down to mps/my = 0.65. We check the effect of extrapolation by performing a fit of A^t = 4 data using only 
the points in the range mps/my = 0.725-0.972. We find that the difference between this fit and our full fit is less 
than 1% for mps/my = 0.65-0.7. 

Finally, we normalize the temperature T/my by the pseudo-critical temperature Tpc/my on the same line of 
constant physics. Results of the procedure above for the lines of constant temperature are shown by dashed lines in 
Fig. for T/Tpc = 1.0, 1.2, 1.4, 1.6, 1.8 and 2.0, where KtiNt = 4) is used to set Tpc. We observe that the Kt line for 
A't = 6 is slightly discrepant from the T/Tpc — 1.5 line; this deviation represents scaling violation in Tpc/my. 



V. BETA FUNCTIONS 



The renormalization group flow along the lines of constant physics is described by the beta functions. Their precise 
values are required in a calculation of the energy density e/T^ discussed in Sec. VL In this section, we study the beta 
functions 



5/3 



da 



dK 
da 



(10) 

for fixed values of mps/my, using results for mpga and my a at zero temperature. 

Since my is a physical quantity which we can take as independent of the lattice spacing a, the derivatives 
and a^f- with fixed mpg/my can be replaced by mya ^^^^^^ and myO g^-^^^^ . Naively, these quantities may be 
determined in the following way. First, one fits the values of mpsa and my a measured at each {(3,K) by a suitable 
fit function, and differentiate the function in terms of B and K. The derivatives at^^ ^ and ^ can be calculated 

by solving 



d0 dK \ / a(mva) a(mps/mv) 

9(mva) a(mva) \ _ / 3/3 

9/3 dK I ^ I a(mva) a(mps/mv) 

d(m-ps/mv) 9(mps/mv) / \ dK dK 



(11) 



In practice, however, we find that there exists a region where the matrix on the right hand side becomes almost 
singular, so that the inverse cannot be solved reliably. In particular, when quarks are heavy, because mps/my is 
always close to one, its derivatives in terms of /3 and K cannot be determined precisely. 

This leads us to adopt the following alternative method. We determine myC g^^^ and myO g^^^ directly from 
the inverse functions /3(mya, mps/wy) and K{mya,mps/my). In Fig. ||, we plot results for (mya, mpg/my) at 
zero temperature. To determine and at a point, say, ((mya)o, ('tips/tov)o), we fit data for (/3, if), i.e., 
the values specifying the simulation points, by a power function expanded in terms of (A(TOva), A(r7ips/my)) = 
{mya — (mya)o, (mps/my) — imps/my)o). We employ the following general fit ansatz up to the third power, 

/3 = c/30 + Ci3i{Amya} + Cf32{Amya}'^ + C/33{Amya}^ + C/34{A{mps/my)} 

+C05{A{mps/my)}{Amya} + Cf36{A{mps/my)}{Amyay + Cf3j{A{mps/my)y 
+Ci38{A{mps/my)y {mya - (mya)o} + C/39{A(mps/r7iv)}^ (12) 



5 



K = ckq + CKi{AmYa} + CK2{Ar7iva}^ + c/f3{Amva}^ + CK4{A{mps/my)} 

+c_ft:5{A(mps/mv)}{ATOpsa} + c_ft:6{A(mps/mv)}{Amva}^ + CK7{A(mps/TOv)}^ 
+Ci<-8{A(mps/mv)}^{Amva} + CRgiM'^Ps/mv)}^ . (13) 

The fit range is determined for each {{mya)o, {m-ps/'n^Y)o) separately: The fit range in mps/my is fixed such that 
X^/-^df is minimized, under the condition that the number of fitted data is larger than 30 to avoid artifacts from 
statistical fluctuations. For the fit range in mya, we inchide aU data except for the points mya < 2.3 or > 0.11, 
which are far from the region we study. 

From the fit results, we calculate the beta functions by differentiating /3 and K in terms of mya, with mps/my 
fixed, 

^vaTT, r = myacpi, (14) 

o[mya) 

dK , , 

mya-— :- = myacKi- (15) 

a(mya) 

The results are plotted in Figs. ^ and |o|. In Fig. ^, the one-loop perturbative value of "ivaa(^a) m^LSsless quark 
is shown by a solid line near the right edge of the plot. Our results appear to gradually approach this value in the 
large (3 limit. We also see that mvQ g(-^^^-| for small mps/my approaches zero at large (3, as we expect from the fact 
that A'c ^ 1/8 as /3 — > oo. 

Another application of the fits (|l2|) and (|l^) is the determination of the lines of constant physics and the temperature 
measured by the vector meson mass, T/my = 1/ (Ntmya), discussed in the previous section. 

VI. EQUATION OF STATE 

The energy density and pressure are defined by 



1 d\nZ d\nZ 



very. 

m 



where Z, T and V are the partition function, temperature and spatial volume, respectively. We calculate these 
quantities as a function of temperature along the lines of constant physics obtained in Sec. 

A. Integral method in full QCD 



We compute the pressure by the integral method This method is based on the formula p = — /, with 

/ = {—T/V) In Z the free energy density, valid for large homogeneous systems. The pressure is then given by 



sub 



where d^ — (d/3', dK') is the line element in the (/?, K) plane, and (• • is the expectation value at finite temperature 
with the zero-temperature value subtracted. The starting point of the integration path should be chosen in the low 
temperature phase where p ~ 0. In actual simulations, for setting the starting point of the integration path, we 
quadratically extrapolate results for the integrand near the low temperature phase to zero. 
For our action (|^) and (||), the derivatives in ( [l7|) are given by 



dlnZ I dS 



dp \ 8(3 



NfNt 



dlnZ _ I dS\ 
dK ~ \dK/ 



NfNfN, 



-^/%^^( y]'^^'''''o^,F^,{x)D-\x,x))\ (18) 
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+ {l + l^.)Ul{x)D-\x,x + fL)])-cswi '^^''"'"^^^.uF^Ax)D-\x,x)\Y (19) 

where Ng is the spatial lattice size and Nf — 2 denotes the number of flavors. We evaluate the quark contributions, 
1^ and using the noisy source method. In order to select the type of noise, we have tested Z{2), U{1) and a 

complex Gaussian noise with a test run on an 8^ x 4 lattice. We find that the U{\) and Z{2) noises show faster 
convergence than the Gaussian noise in the number of noise ensembles. The difference between the U{1) and Z{2) 
cases is small, while the U{\) noise shows slightly faster convergence in this test. From this result, we have adopted 
the U(\) noise in this study. 

The integral method was originally developed for a pure gauge system, for which the parameter space is one- 
dimensional. In our case, the parameter space is two-dimensional. Therefore, the integration path for the pressure 
( p7| ) is not unique. We have performed a series of test runs on an 8^ x 4 lattice, and have confirmed the integration 
path independence . Details of the test are presented in Appendix ^ We shall present results of a similar test in 
our production runs below. 

We also find from this test that the integration paths in the K direction (constant 0) lead to smaller errors for the 
final values of the pressure than those in the (3 direction (constant K). Therefore, we choose the integration paths in 
the K direction starting from the region of small values of K and moving towards the chiral limit. Our simulation 
points on the x 4 and x 6 lattices are shown by "-f" and "x" in Fig. ||. 



B. Equation of state for Nt 



In Fig. |ri|, we show the results for the pressure derivative 

d{plT^) _ .1 1 dS\ 



dK 



sub 



obtained on an A^t = 4 lattice. Measurements are made with five noise ensembles at every trajectory. The bin size for 
the jack-knife errors is set to 10 trajectories from a study of bin size dependence. Numerical values for the derivative 
are summarized in Table [v^ . Interpolating the data by a cubic spline method, we integrate in the K direction to 
obtain the pressure presented in Fig. |l2|. 

We also compute the derivative in the /3 direction, 

d(p/n-^N^/^^\ (21) 



as shown in Fig. n3| and Table VII. We observe that the results for this derivativ e are noisier than those for the K 



derivative in Fig. |ll| . This is the underlying reason for the fact commented in Sec. VIA that the integral paths in the 
K direction lead to smaller errors in the pressure. 

In Fig. |l^, we replot the pressure from the K integration as a function of /3, and compare it with the slope ( [2l| ) 
computed independently from the ^ data. The latter data for the slope, shown by short lines, are tangential to the 
pressure curve, confirming the integration path independence of results for the pressure. 

The pressure data shown in Fig. |l^ or [l4| arc not quite useful yet. We wish to compute the pressure on a line 
of constant physics as a function of temperature normalized by the pseudo-critical temperature Tpc on the same 
line of constant physics. The necessary change of parameters from {(3,K) to (mps/my, T/Tpc) is achieved with the 
interpolations performed in Sees. IV and 0. 



The pressure p/T^ as a function of T/Tpe is given in Fig. |T| for mps/my = 0.65, 0.7, 0.75, 0.8, 0.85, 0.9 and 0.95. 
In this figure, symbols denote the values on the integration path along the K direction &t (3 — 1.80, 1.85, 1.90, 1.95, 
2.00, 2.10 and 2.20, i.e., for each /3, the pressure in Fig. O at the values of K corresponding to the given values of 



mps/rnv- The values of T/Tpc at those points are determined from mya, as discussed in Sec. IV. To interpolate 

as 



these symbols in the direction of (3 {i.e., of T/Tpc ), we use the results for the slopes ^ shown in Fig. 14 and adopt 



a cubic ansatz. 

We observe in Fig. |l^ that the pressure depends only weakly on the quark mass once the ratio mps/my falls below 
« 0.8. In the heavy quark limit mps/my = 1, the pressure should coincide with the pure gauge value on a lattice 
with the same size, which is denoted by a dashed line While the pressure decreases for mps/my = 0.80-0.95, the 
values at mps/my — 0.95 are still large compared to those of the pure gauge system for A''^ = 4. 

In Fig. |l^, we also mark the high-temperature Stefan-Boltzmann (SB) values by solid horizontal lines, both for 
Nt = 4: and in the continuum. The lattice value is evaluated from the free energy density in the SB limit, to be in 
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parallel with the integral method adopted in numerical simulations. Some details of this computation are described in 
Appendix |b[ We observe that the pressure overshoots the SB value in the continuum limit, and appears to gradually 
increase toward the SB value for the TVt = 4 lattice at high temperatures. Another point to note is that the large SB 
value on an Nt — i lattice [||j3^ is dominated by the quark contribution. While the integral method does not allow a 
separate evaluation of the two contributions, a comparison of the two-flavor result and that of the pure giuon theory 
1^ (dashed line) shows that the situation should be similar at finite temperatures. 
To compute the energy density, we use the following expression for the interaction measure e — 3p: 



-3p J 1 ds\ 



sub 

^ r Pia I 1 



d(5 I 1 dS\ dK I 1 dS \ 



(22) 



Applying the beta functions calculated in Sec. ^ we find the results for (e — 3p)/T^ shown in Fig. |l6|. The meaning 
of symbols is the same as in Fig. 

Combining Figs. |l| and |l| for p/T^ and (e - 3p) /T"^, we obtain the energy density plotted in Fig This quantity 
also overshoots the SB value in the continuum limit. In contrast to the case of pressure, the energy density in the 
high temperature phase is quite constant as a function of temperature. 

Our results for pressure and energy density allow us to calculate the speed of sound defined by 

^ |. (23) 



e 



We compute the derivative from a quadratic fit of p as a function of e using three data points along the lines of 
constant physics. The results for iVj = 4 are plotted in Fig. |l^ where errors are estimated by error propagation from 
those of p and e. We omit results at small temperatures with T/Tps < 0.9, since there the magnitude of e is small 
in comparison with its error. The speed of sound rapidly increases just above the critical temperature, and almost 
saturates the SB value when T/Tpc > 1.5. 



C. Equation of state for Nt = 6 

The simulation points for Nt — 6 lattices at = 1.95, 2.00, 2.10, 2.20, and 2.30 are marked by crosses in Fig. ||. 
Raw data for the derivatives d{p/T'^) / dK and d{p/T'^) / d(3 are shown in Fig. |l^ and Fig. |2^, respectively. Since the 
statistical errors for Nt = 6 results are larger than those for Nt — 4, the spline interpolation does not work well for 
Nt — 6. Therefore, we interpolate the pressure derivatives by straight lines. The rest of data analyses parallel those 
for the case of Nt = 4. 

The pressure for A't = 6 on the lines of constant physics are plotted in Fig. ^ by open symbols as a function of 
T/Tpc, together with the results for Nt = 4 (filled symbols). Figure ^ shows the interaction measure (e — 3p)/T^ 
as a function of T/Tpc, calculated from results in Figs. |l^ and pO| together with the beta functions of Figs. p| and [ic| . 
Combining the results for pressure and interaction measure, we obtain the energy density presented in Fig. 12^. 

These figures show that both the pressure and energy density decreases as the lattice spacing becomes smaller from 
A't = 4 to 6, and the values at high temperatures become closer to the continuum SB limit. On the other hand, both 
at At = 4 and 6, the energy density is smaller than the SB values for the corresponding A't, and an approach to the 
lattice SB value toward high temperatures is not apparent in our data. A similar deviation of energy density from the 
SB value at finite At has been reported in Ref. |^ for the case of the SU(3) pure gauge theory using the RG improved 
gauge action (^). Further study is necessary to examine if deviations remain toward the limit of high temperatures. 

We also observe for both the pressure and the energy density that the dependence on the quark mass is quite small 
for mps/my « 0.65-0.8. A weak quark mass dependence appears only at mps/my >0.9 for A't = 4 (the errors are 
still too large to conclude a quark mass dependence for the Nt = 6 data). This result may not be surprising since 
hadron mass results in our zero-temperature simulations ||28|| show that the renormalized quark mass at /z = 2 GeV 
in the MS scheme at mp^/ra\ w 0.65-0.8 equals m^^^{fj, — 2GeV) w 100-200 MeV, which is similar in magnitude to 
the critical temperature Tc ~ 170 MeV estimated for two-flavor QCD |jl|. For comparison, finite mass corrections for 
free fermion gas only amount to 7% when the temperature equals the fermion mass mf, and exceed 50% only when 
TO//T>3. 

In a previous study using the standard staggered quark action, it was reported that the energy density e/T^ for 
A't = 4 overshoot the SB value forming a sharp peak just above Tc, while the results for Nt = 6 show no peak 
With the improved Wilson quark action, we do not observe an overshoot of the energy density both at A't = 4 and 
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6. Similar absence of the peak of energy density is reported also from an improved staggered quark action at iVt = 4 
when a contribution proportional to the bare quark mass is removed We think it likely that the overshoot 

observed with staggered quark action at A't = 4 is a lattice artifact. With the staggered quark action, the energy 
densities for rUq/T — 0.075 and 0.15 at iVt = 6 are found to be consistent with each other within the errors [||. This 
result is similar to our finding of small quark mass dependence in the EOS. 

Our present data for Nt — 4: and 6 show a 50% decrease both in the pressure and energy density, which is too large to 
attempt a continuum extrapolation. On Nt — 6 lattices, however, the magnitude and temperature dependence of the 
two quantities are quite similar between our improved Wilson quark action and the staggered quark action. Together 
with the fact that the Nt = 6 results are close to the continuum SB limit at high temperatures, the approximate 
agreement of EOS between two different types of actions may be suggesting that the Nt = 6 results are not far from 
the continuum limit. This expectation is also supported by the Nt dependence of the SB value on the lattice. The 
value for Nt — 6, which is 50% too large compared to the continuum limit, reduces by 30% so that for Nt = 8 the 
lattice SB value is within 20% of the continuum limit. Thus we expect that a precise continuum extrapolation will 
be possible if additional data points at A^t = 8 are generated, which we leave for future work. 



VII. CONCLUSION 



We have presented first results for the equation of state in QCD with two flavors of dynamical quarks using a 
Wilson-type quark action. In order to suppress large lattice artifacts observed with the standard Wilson quark action, 
we have adopted a clover-improved form of the action and an RG- improved gluon action. Two temporal lattice sizes, 
TVt = 4 and 6, are studied to examine the magnitude of finite lattice spacing errors. 

We have calculated the energy density and the pressure as functions of temperature along the lines of constant 
physics, which are identified through the mass ratio mps/jTiv- As a part of the analysis to work out these lines, we 
have also computed the beta functions in the parameter space K). 

We found that the quark mass dependence of EOS is small over the range mps/rriY « 0.65-0.8. While the physical 
point mps/wv = 0.18 is still far away, the observed independence on the quark mass suggests that our result for the 
EOS is close to those at the physical point except in the vicinity of the chiral transition point where a singular limit 
according to the 0(4) critical exponents is expected. 

Our results for the pseudo-critical temperature in unit of the vector meson mass, Tpc/rriY, show an agreement within 
about 10% between the temporal lattice sizes A'^f = 4 and 6. On the other hand, the pressure and the energy density 
decreases substantially, showing the presence of large scaling violation in the results for Nt = 4. An encouraging 
indication, however, is that results on the Nt = 6 lattice are close to the continuum Stefan-Boltzmann value at high 
temperatures. We also note that the values and the temperature dependence of EOS at iVt = 6 are quite similar to the 
previous results from staggered quark action [||. These may be suggestions of a possibility that precise calculations 
of EOS are realized on lattices with temporal sizes Nt not much larger than 6. 
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APPENDIX A: A TEST OF THE INTEGRAL METHOD IN FULL QCD 



We compute EOS by the integral method [g7| described in Sec. VIA. In order to test the method, we perform a 
series of test runs to calculate the pressure on an 8^^ x 4 lattice. For subtraction of the zero temperature part, we also 
measure the same operators on an 8"^ lattice. Simulation points are shown by stars in Fig. We generate 500 HMC 
trajectories at each point. 

We first check the influence of the spatial volume on EOS. In order to avoid systematic errors from numerical inter- 
polation and extrapolation needed for numerical integrations, we first compare the values of the integrand d{p/T^)/dK 
for Ns/Nt = 2 (8^ X 4 lattice) with those for Ng/Nt = 4 from our main runs on a 16'^ x 4 lattice. From Fig. [23, we 
find that the diffirence between Ng/Nt = 2 and 4 is comparable with statistical fluctuations. Integrating out these 
values, we obtain Fig. We find that the two results agree well with each other — the slight discrepancy of the 
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P = 2.2 data at small K seems to be caused by a longer extrapolation to the point d{p/T^) / dK = for the stating 
point of the integration. 

We then test the integration path independence in the integral method. We study three integration paths in the 
parameter space {f3,K), shown in Fig. g^. At the crossing points, the results for the pressure from different paths 



should coincide with each other. The results for p/T"^ obtained from these paths are summarized in Fig. 27, The left 
figure is obtained by integrating in the /3 direction at K ^ 0.13, while the right figure is computed from the integration 
paths in the K direction at (3 = 2.1 and 2.2. 

We find that p/T** at K) = (2.1, 0.13) and (2.2, 0.13) in the two figures agree well with each other. This confirms 
the path-independence of the pressure. 

We also note that the paths in the K direction lead to much smaller errors in the pressure than the path in the /3 
direction. We therefore adopt paths in the K direction in the production runs discussed in the main text. 

APPENDIX B: STEFAN-BOLTZMANN LIMIT OF PRESSURE BY THE INTEGRAL METHOD 



In the calculation of pressure discussed in Sec. VI, we employ the integral method, in which the negative of free 
energy density — / = {T/V)\nZ is identified with the pressure. In order to compute the Stefan-Boltzman value to be 
compared with the pressure from the integral method, we should compute the free energy density in the free gas limit. 
In this appendix, we describe our calculation of the free energy density for the case of our improved lattice action. 

To calculate the partition function Z, we expand the link variable as 

U^ix) = exp{igA^{x)} = exp{igA''^{x)Ta} (Bl) 
and perform a Fourier transformation 

A^,(x) = ^^=y e*'=(^+'^/2)A;;(fc), (B2) 



where 



= j^=0,±l,---,iV,/2 forAi = l,2,3 (B3) 

k^. = = 0, ±1, • • • , Nt/2 for ^i = 4. (B4) 



Fixing the gauge to the lattice Lorentz gauge by adding the gauge fixing term 



S,f = -J2Tr(j2^Mx)] , 



(B5) 



the free part of the gauge action (||) is given by 
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where = 2sin(fc^/2), (fc) = i(fcpA«(fc) - KAl{k)), and 



q^. = 1 - ciCkl + kl); iorii^v (B7) 
q,^u = 0; for /i u. (B8) 



The free part of the ghost term corresponding to the gauge fixing (pa) is given by 



4h = E KKt(~k)7^''(k), (B9) 



where rj and rj are the ghost fields. 

The partition function for the gauge part can be calculated as 
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^A&lD]l^{k)Dl^^(k) X (const.) 



with 



p 



Consequently, we obtain the gauge part of the unnormahzed free energy density 



1 



■ hi Z„ = 



Y.' 



5 ln(detD;„'(t))- In 



where means a sum except for the zero mode. 
The free part of the quark action is given by 



S, = -NfY.Y.^ik)D,q^ik), 



a=l k 

4 



Dq{k) = 2Kj2^osikp) - l-2iif^7^sin(/c^), 

fj.— l M— 1 



where 



jV,-0,±l,---,iV,/2 iov^l = 1,2,3 



2%tlZH), j^ = 0,±l,...,7V,/2 forM = 4. 



The partition function for the quark part is obtained as 



Z^= (^dct D,{k)^ 



3Nf 



(const.). 



(BIO) 

(Bll) 
(B12) 

(B13) 



(B14) 
(B15) 



(B16) 
(B17) 



(B18) 



det Dqik) 



1-8K + AK^sin^{kp/2)\ + AK^J^^^^^^i^f 



We then obtain the quark part of the unnormahzed free energy density [p2[ 



NfNt 



(B19) 



Numerical calculations of the normalized energy density and pressure are performed using the equations, 



pa 



(B20) 



(B21) 
(B22) 



where /t^h and J^Iq are the free energy at zero temperature calculated on an lattice. At the right edge of the 
figures [l5|^p^, pl| and |2^, we show the results for the cases of 16^ x 4 and 16'^ x 6 lattices, in the massless quark limit, 
K = l/8. 



fig) 
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p 


K 


traj. 


therm. 


1.80 


0.1300-0.1450 


500-2000 


200-500 


1.85 


0.1250-0.1440 


500-1900 


200-300 


1.90 


0.1250-0.1425 


500-2000 


200-400 


1.95 


0.1200-0.1410 


500 2000 


200 


2.00 


0.1150-0.1390 


500-2000 


200-300 


2.10 


0.0900-0.1375 


500-1000 


200-900 


2.20 


0.0700-0.1365 


500 


200 


TABLE I. Simulation parameters on a 16"^ x 4 lattice. 


P 


K 


traj. 


therm. 


1.95 


0.1350-0.1410 


1000 


200 


2.00 


0.1300-0.1395 


800-1500 


200 


2.10 


0.1200-0.1375 


500-1000 


200-300 


2.20 


0.1100-0.1365 


500-1000 


200-300 


2.30 


0.1000 O.l.'ifiO 


m) vm 


200 250 




TABLE II. Simulation parameters on a 16 


^ X 6 lattice. 
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K 


traj. 


therm. 


1.80 


0.1300-0.1450 


200 


200-500 


1.85 


0.1250-0.1440 


200-300 


100-300 


1.90 


0.1250-0.1425 


200 


200 400 


1.95 


0.1200 0.1410 


200-300 


100 400 


2.00 


0.1150 0.1395 


200-300 


100-200 


2.10 


0.0900-0.1375 


300 


200-550 


2.20 


0.0700-0.1365 


200-300 


100-200 


2.30 


0.1000-0.1360 


200 


100 



TABLE III. Simulation parameters on a 16* lattice. 
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0.1250 


1.2666(25) 


1.3323(29) 


0.9507(9) 


2.20 


0.1275 


1.1077(22) 


1.1788(26) 


0.9397(11) 


2.20 


0.1300 


0.9289(26) 


1.0167(38) 


0.9137(21) 
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2.20 


0.1325 




0.7456(33) 


0.8442(34) 


0.8831(29) 


2.20 


0.1350 




0.5014(31) 


0.6307(55) 


0.7950(51) 


2.30 


0.1100 




1.9371(18) 


1.9655(20) 


0.9855(2) 


2.30 


0.1150 




1.6957(20) 


1.7314(22) 


0.9794(4) 


2.30 


0.1175 




1.5721(20) 


1.6102(23) 


0.9763(4) 


2.30 


0.1200 




1.4292(20) 


1.4692(23) 


0.9728(5) 


2.30 


0.1225 




1.2959(20) 


1.3397(19) 


0.9673(5) 


2.30 


0.1250 




1.1539(29) 


1.2074(29) 


0.9556(11) 


z.oU 


U.lz /o 










2.30 


0.1300 




0.8117(25) 


0.8812(30) 


0.9211(13) 


2.30 


0.1325 




0.6297(29) 


0.7110(47) 


0.8857(29) 


2.30 


0.1340 




0.4867(33) 


0.5932(40) 


0.8205(51) 


2.30 


0.1350 




0.3807(36) 


0.5080(53) 


0.7493(69) 




TABLE IV. Results for mpso, 


mvo and mps/mv on a zero temperature lattice of size 16*. 
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Kt{Nt=A) 


mps /mv 




Kt{Nt = 6) mps/mv 


Tpc/my 


1.600 


0.1543(10) 


0.346(153) 


0.217(11) 






1.650 


0.1533(10) 










1.700 


0.1510(10) 


0.396(170) 


0.234(17) 






1.800 


0.1445(14) 


0.690(92) 


0.211(15) 






1.850 


0.14019(18) 


0.7905(60) 


0.1917(20) 






1.900 


0.13621(15) 


0.8525(39) 


0.1801(12) 






1.925 


0.13417(23) 










1.950 


0.13040(97) 


0.9051(64) 


0.1572(62) 






2.000 


0.12371(73) 


0.9450(36) 


0.1398(29) 


0.13861(21) 0.725(16) 


0.2086(53) 


2.100 


0.10921(43) 


0.9790(13) 


0.1114(09) 


0.13365(40) 0.8635(78) 


0.1753(58) 


2.200 








0.12539(25) 0.9481(19) 


0.1275(15) 


2.300 








0.11963(15) 0.9724(12) 


0.11145(62) 



TABLE V. Finite temperature transition/crossover point Kt for Nt = 4 and 6. Results for mps(T = 0)/mv{T = 0) and 
Tpc/my {T = 0) interpolated to the Kt point are also listed. 
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1^ ^ ~Nf7h(^) 







16^ X 4 


16^ X 6 


16* 


1.80 


0.1300 


-5.3622(47) 




-5.4232(59) 


1.80 


0.1350 


-4.8166(79) 




-4.9341(51) 


1.80 


0.1375 


-4.3207(63) 




-4.5143(54) 


1.80 


0.1400 


-3.5759(65) 




-3.9051(93) 


1.80 


0.1425 


-2.2016(162) 




-2.9788(112) 


1.80 


0.1440 


-0.2992(227) 




-2.1623(102) 


1.80 


0.1450 


2.1328(221) 




-1.4493(102) 


1.85 


0.1250 


-4.7935(43) 




-4.8429(33) 


1.85 


0.1300 


-4.3596(47) 




-4.4437(33) 


1.85 


0.1350 


-3.4384(56) 




-3.6678(50) 


1.85 


0.1375 


-2.6171(80) 




-3.0360(71) 


1.85 


0.1400 


-0.8984(161) 




-2.0797(142) 


1.85 


0.1425 


2.6710(146) 




-0.6809(104) 


1.85 


0.1440 


4.5010(140) 




0.7871(132) 


1.90 


0.1250 


-3.9170(71) 




-3.9860(32) 


1.90 


0.1300 


-3.2151(80) 




-3.3875(53) 


1.90 


0.1325 


-2.6335(101) 




-2.9136(34) 


1.90 


0.1350 


-1.5613(124) 




-2.2669(69) 


1.90 


0.1375 


0.3121(130) 




-1.3867(70) 


1.90 


0.1400 


2.6656(93) 




-0.0358(137) 


1.90 


0.1425 


4.9308(89) 




2.0916(163) 


1.95 


0.1200 


-3.5434(50) 




-3.6066(38) 


1.95 


0.1250 


-2.9842(63) 




-3.1191(60) 


1.95 


0.1275 


-2.5015(123) 




-2.7719(46) 


1.95 


0.1300 


-1.7833(65) 




-2.3076(73) 


1.95 


0.1325 


-0.6182(99) 




-1.6726(52) 


1.95 


0.1350 


0.8561(76) 


-0.8309(61) 


-0.8373(61) 


1.95 


0.1375 


2.4815(64) 


0.3734(92) 


0.3009(84) 


1.95 


0.1390 




1.3834(87) 


1.2217(86) 


1.95 


0.1400 


4.2807(63) 


2.3321(134) 


2.0463(106) 


1.95 


0.1410 


5.0293(99) 


3.6109(165) 


2.9673(135) 


2.00 


0.1150 


-3.2222(42) 




-3.2710(27) 


2.00 


0.1200 


-2.7735(49) 




-2.8989(40) 


2.00 


0.1250 


-1.8085(56) 




-2.2919(32) 


2.00 


0.1275 


-1.0607(82) 




-1.8253(43) 


2.00 


0.1300 


-0.0836(73) 


-1.2623(54) 


-1.2600(49) 


2.00 


0.1325 


0.9962(91) 


-0.5093(40) 


-0.5103(75) 


2.00 


0.1350 


2.2580(69) 


0.5282(49) 


0.4820(82) 


2.00 


0.1375 


3.7258(78) 


2.0073(86) 


1.8917(77) 


2.00 


0.1385 




2.9849(14) 


2.5623(111) 


2.00 


0.1390 


4.6687(79) 


3.4516(133) 


2.9676(112) 


2.00 


0.1395 




4.0484(106) 


3.3763(124) 


2.10 


0.0900 


-2.9462(25) 




-2.9600(17) 


2.10 


0.1000 


-2.8643(27) 




-2.8993(14) 


2.10 


0.1050 


-2.7054(31) 




-2.7710(25) 


2.10 


0.1100 


-2.3368(50) 




-2.5538(26) 


2.10 


0.1150 


-1.8240(60) 




-2.1962(21) 


2.10 


0.1200 


-1.0100(60) 


-1.6474(34) 


-1.6509(29) 


2.10 


0.1250 


0.0614(57) 


-0.8130(62) 


-0.8219(32) 


2.10 


0.1300 


1.6577(59) 


0.5028(66) 


0.4794(45) 


2.10 


0.1325 




1.5077(73) 


1.3864(71) 


2.10 


0.1340 




2.3734(77) 


2.0926(79) 


2.10 


0.1350 


3.9084(51) 


2.9291(67) 


2.5742(52) 


2.10 


0.1365 




3.9064(99) 


3.4686(71) 


2.10 


0.1375 


5.3108(39) 


4.6615(65) 


4.1983(60) 


2.20 


0.0700 


-2.3029(21) 




-2.3022(16) 


2.20 


0.0800 


-2.4132(27) 




-2.4469(16) 
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2.20 


0.0900 


-2.3649(31) 




-2.4697(21) 


2.20 


0.1000 


-2.0985(35) 




-2.3012(25) 


2.20 


0.1100 


-1.4253(42) 


-1.8043(21) 


-1.8065(20) 


2.20 


0.1200 


-0.0533(45) 


-0.6958(26) 


-0.6954(24) 


2.20 


0.1225 




-0.2433(30) 


-0.2560(41) 


2.20 


0.1250 


1.1380(43) 


0.3244(35) 


0.2736(27) 


2.20 


0.1275 




1.0372(42) 


0.9297(37) 


2.20 


0.1300 


2.7957(50) 


1.9547(36) 


1.7599(28) 


2.20 


0.1325 


3.8577(40) 


3.0423(41) 


2.7646(49) 


2.20 


0.1350 


5.0923(39) 


4.4078(41) 


4.0746(41) 


2.20 


0.1365 


5.8931(36) 


5.3480(43) 


5.0113(52) 


2.30 


0.1000 




-1.8180(16) 


-1.8210(21) 


2.30 


0.1100 




-1.2098(17) 


-1.2178(21) 


2.30 


0.1150 




-0.6777(21) 


-0.6944(21) 


2.30 


0.1175 




-0.3082(20) 


-0.3570(29) 


2.30 


0.1200 




0.1051(23) 


0.0439(28) 


2.30 


0.1225 




0.6267(24) 


0.5338(29) 


2.30 


0.1250 




1.2434(23) 


1.1185(25) 


2.30 


0.1275 




2.0087(33) 


1.8252(21) 


2.30 


0.1300 




2.9362(34) 


2.7110(40) 


2.30 


0.1325 




4 0641 C^^) 


3 8n77('3S'l 


2.30 


0.1340 




4.8553(29) 


4.6173(29) 


2.30 


0.1350 




5.4316(35) 


5.1741(50) 


2.30 


0.1355 




5.7271(42) 


5.4944(29) 


2.30 


0.1360 




6.0362(32) 


5.8475(45) 



TABLE VI. Derivative (H) for 16^ x 4, 16* x 6 and 16* lattices. 
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2.20 


0.0900 


11.8419(15) 




11.8180(10) 


2.20 
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2.20 
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0.1150 
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12.3021(9) 


2.30 


0.1200 




12.3234(6) 


12.3187(9) 


2.30 


0.1225 




12.3420(6) 


12.3348(8) 


2.30 


0.1250 




12.3604(7) 


12.3550(11) 


2.30 


0.1275 




12.3812(8) 


12.3722(7) 


2.30 


0.1300 




12.4034(9) 


12.3947(9) 


2.30 


0.1325 




1 2 4241 (H) 


12 4^HO(^o) 

iZj.^_L(JWl -L\J J 


2.30 


0.1340 




12.4368(8) 


12.4337(8) 


2.30 


0.1350 




12.4464(8) 


12.4401(8) 


2.30 


0.1355 




12.4497(9) 


12.4454(8) 


2.30 


0.1360 




12.4550(7) 


12.4509(11) 



TABLE VII. Derivative (If ) for x 4, 16^ x 6 and 16* lattices. 
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FIG. 1. Pseudoscalar meson mass squared as a function of 1/K. Filled symbols are obtained on a 16* lattice. Open symbols 
are from Ref. p^ . 




1/K 

FIG. 2. Same as Fig. hi, but for vector meson mass. 
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FIG. 3. Phase diagram and simulation points on 16'^ x 4, 16^ x 6 and 16* lattices. 
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FIG. 4. Polyakov loop obtained on a 16 x 6 lattice. 
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FIG. 5. Polyakov loop susceptibility obtained on a 16^ x 6 lattice. 
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FIG. 6. Pseudocritical temperature in units of vector meson mass as a function of (mps/mv)^. Dashed lines are interpola- 
tions, separately, for iVt = 4 and 6 data. 
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FIG. 7. Lines of constant physics and of constant temperature. Solid lines are mps/mv constant lines, and dashed lines are 
T/Tpc constant lines for Nt — 4. The values of T/Tpc for the dashed lines are given on the right edge of the figure. 
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FIG. 8. mps/mv versus mya obtained on T = lattices. Solid lies are (3 constant lines, and dashed lines are K constant 
lines. Meaning of symbols are the same as those in Fig. |l| 
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FIG. 9. mya 



on nivs/'mY constant lines. One- loop result from perturbation theory at rUq = is also shown. 
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FIG. 10. 
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on mps /mv constant lines. 
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FIG. 13. Derivative of pressure with respect to /3 as a function of on a 16^ x 4 lattice. 
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FIG. 14. Pressure as a function of /? on a 16^ x 4 lattice. Short Unes denote the slope of the pressure for /3 direction. 
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FIG. 15. Pressure on a 16"^ x 4 lattice as a function of T /Tpc. The dashed curve shows the pressure for pure gauge theory 
with the RG-improved action on a 16^ x 4 lattice Q. 
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FIG. 16. (e - Sp)/?* on a 16^ x 4 lattice as a function of T /Tpc- 
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FIG. 17. Energy density on a 16'^ x 4 lattice as a lunction of T/Tpc- The dashed curve shows the energy density for pure 
^auge theory with the RG-improved action on a 16^ x 4 lattice Q. 
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FIG. 18. Speed of sound squared on a 16^ x 4 lattice as a function of T/Tpc. 
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FIG. 19. Derivative of pressure with respect to as a function of if on a 16^ x 6 lattice. 
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FIG. 20. Derivative of pressure with respect to /3 as a function of X on a 16^ x 6 lattice. 
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FIG. 21. Pressure on 16^ x 4 (filled symbols) and 16^ x 6 (open symbols) lattices as a function of T/Tj 




FIG. 22. (e - 3p)/r* on 16^ x 4 (filled symbols) and 16^ x 6 (open symbols) lattices as a function of T/Tj 
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FIG. 24. Integration paths for test simulations on an 8^ x 4 lattice. Stars represent the simulation points. 
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FIG. 25. d{p/T'^)ldK on 8^ x 4 and 16^ x 4 lattices. 
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FIG. 26. The pressure p/T* on 8^ x 4 and 16^ x 4 lattices. 
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